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We study the effects of size polydispersity on the gas-Hquid phase behaviour of 
mixtures of sticky hard spheres. To achieve this, the system of coupled quadratic 
equations for the contact values of the partial cavity functions of the Percus-Yevick 
solution is solved within a perturbation expansion in the polydispersity, i.e. the nor- 
malized width of the size distribution. This allows us to make predictions for various 
thermodynamic quantities which can be tested against numerical simulations and 
experiments. In particular, we determine the leading-order effects of size polydis- 
persity on the cloud curve delimiting the region of two-phase coexistence and on the 
associated shadow curve; we also study the extent of size fractionation between the 
coexisting phases. Different choices for the size-dependence of the adhesion strengths 
are examined carefully; the Asakura-Oosawa model of a mixture of polydisperse col- 
loids and small polymers is studied as a specific example. 
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I. INTRODUCTION 

In the context of soft matter, a number of systems are known to display a combination 
of a very steep repulsion and a short range attraction. This includes for instance polymer- 
coated colloidsi*^, globular proteins^ and microemulsions^. In spite of the notable differences 
in the details of the interactions among these systems, most of the common essential features 
are captured by a paradigmatic model known as the adhesive or sticky hard sphere model. 
Sticky hard spheres are impenetrable particles of diameters {cTj} with an adhesive surface. 
The simplest way of describing the adhesion properties, in the framework of atomic fluids, 
was originally proposed by Baxter— in terms of a potential where energy and length scales 
were combined into a single parameter, thus defining the so-called sticky hard sphere (SHS) 
potential. Baxter showed that for this model the Ornstein-Zernicke integral equation de- 
termining the correlation functions in the liquid state admitted an analytic solution within 
the Percus-Yevick (PY) approximation. Together with his collaborators, he predicted from 
this solution (via both the compressibility and the energy routes of liquid state theory) 
that the model displays a gas-liquid transition^. This PY solution was soon extended to 
mixtureaSiSiifliii and has since found a number of interesting applications in the area of col- 
loidal suspensionsii^iiSii^iiiii^. When studying the phase behavior of such fluids an important 
issue to deal with is the fact that colloidal particles are generally not identical but may have 
different characteristics (size, charge, chemical species etc). Often, the distribution of the 
relevant parameter is effectively continuous and the fluid is then referred to as polydisperse. 
We will focus in this paper on size polydispersity, i.e. a fluid with a distribution of particle 
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diameters. (A small degree of size polydispersity is in fact required to resolve thermody- 
namic pathologies which occur in the case of strictly equal-sized, i.e. monodisperse, sticky 
hard spheres^S,.) The particle size distribution is fixed when the particles are synthesized. 
Thereafter, only the overall density can be modified by adding or removing solvent, while 
keeping constant all ratios of densities of particles of different size; this traces out a so-called 
"dilution line" in the phase diagram. 

Given the success of the PY closure for the monodisperse SHS model, it is natural to 
try to extend it to the polydisperse case. Unfortunately, the PY approximation is tractable 
only for mixtures of a small number of particle species: the case of a binary mixture can be 
solved analyticallyii, and for mixtures with a limited number of components (10 or fewer) a 
numerical solution is feasible^^. The polydisperse case requires one to keep track of an effec- 
tively infinite number of particles species, one for each size, and cannot be tackled directly. 
An alternative, which we have explored in past work, is to use simpler integral equation 
theories such as the modified Mean Spherical Approximation (mMSA or CO). Between this 
and the Percus Yevick (PY) approximation^i lie a set of increasingly accurate approxima- 
tions, denoted as Cn with n = 1,2,... They are based on a density expansion of the direct 
correlation function outside the hard core and can be shown to improve, order by order, the 
various virial coefficientsi^. These Cn approximations can be extended to the polydisperse 
case with relative ease, provided a particular factorization holds for the matrices appear- 
ing in the solution of Baxter's equations. This has allowed us to perform a comprehensive 
analysis of polydispersity effects on the gas-liquid phase separationi^. 

The tractability of the Cn approximations for the polydisperse SHS model does, however, 
come at the price of lower accuracy. Indeed, for the monodisperse case accurate Monte Carlo 
simulation data recently published by Miller and Frenke l^^i^°i^^ show that the equation of 
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state of the fluid lies very close to the one derived from the energy route of the PY closure. 
Both the CO and CI approximations, on the other hand, yield precise results only within 
a rather limited region of the phase diagram, corresponding to high temperatures or low 
densities^; see Fig. [T] below. 

The above considerations show that another attack on the PY closure for polydisperse 
SHS fluids is worthwhile in order to get accurate predictions for the gas-liquid phase behav- 
ior. Rather than trying to tackle the most general case of a fluid with a potentially wide 
distribution of particle sizes, which for now remains out of reach, we exploit the idea of 
Evans^^ to treat size polydispersity as a perturbation to the monodisperse phase behavior. 
For this method to apply, the size distribution only has to be sufficiently narrow but its 
shape is otherwise arbitrary. Our approach is also of sufficient generality to consider arbi- 
trary dependences of the adhesion strengths on the particle sizes, including those considered 
in previous work on the Cn approximation o^^'^? . Throughout, we consider gas-liquid phase 
coexistence only. It has been argued^i that even in the presence of polydispersity this is 
metastable with respect to phase separation into a colloidal gas and solid. However, the 
latter may be unobservable on realistic timescales when formation of the polydisperse solid 
is hindered by large nucleation barriers2£ or an intervening kinetic glass transition's,; the 
gas-liquid phase splits we calculate will then control the physically observable behaviour. 
Even where the kinetics does allow formation of solid phases, the metastable gas-liquid phase 
behaviour can play a role, e.g. in determining phase ordering pathwaya^i. 

The paper is organized as follows. In section IH] we describe the polydisperse SHS model 
and discuss various routes for predicting the thermodynamics of this system, comparing 
their accuracy for the better understood monodisperse case. In the polydisperse setting one 
needs to model how the strength of the adhesion between two particles depends on their 



size; we discuss some possible choices for this in section ITTTl Section |TVl describes our pertur- 
bation expansion of the PY closure for the weakly polydisperse SHS model. We first define 
the perturbation expansion of the free energy used by Evans (Sec. IIV Ajl and summarize the 
relevant consequences for two-phase coexistence and the attendant size fractionation effects. 
The basic equations that need to be solved in order to determine thermodynamic proper- 
ties within the PY approximation are then described and solved perturbatively fSec. IIVB|) . 
while Sec. IIV(yl derives from this, via the energy route, the excess Helmholtz free energy. 
In section El we evaluate numerically the consequences of polydispersity for two-phase coex- 
istence and fractionation for a number of example scenarios, and compare with alternative 
approximation schemes. Section IVll gives concluding remarks. 



The p-component SHS mixture model is made up of Hard Spheres (HS) of different 
diameters ai, i = 1,2, ... ,p, interacting through a particular pair potential defined via the 
following limit procedure. One starts with a pair interaction potential (pijir) with a hard 
core extending out to distance r = aij = {ai + crj)/2, followed by a square well potential of 
width Rji — cr,,: 



II. THE SHS MODEL 
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Here the dimensionless parameter 
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measures the surface adhesion strength or "stickiness" between particles of species i and j. 
In Eq. (121) the reduced temperature r is an unspecified increasing function of the physical 
temperature T; the coefficients eij specify how stickiness depends on which particle species 
are in contact and are discussed more fully in the next section. The procedure which defines 
the SHS model then consists in taking the "sticky limit" Rij — * aij. The logarithm in 
the initial square well potential is chosen such as to give a simple expression for the 
Boltzmann factor exp[— 0jj(r)], which reduces to a combination of a Heaviside step function 
and a Dirac delta function in the sticky limit. Here and in the following we measure all 
energies in units of k^T, to simplify the notation. 

A fully polydisperse system is obtained from the above discrete mixture by replacing the 
molar fractions Xi = Ni/N, where Ni is the number of particles of species i and the total 
number of particles, with a normalized size distribution function p(cr): 



Here p{a) da is the fraction of spheres with diameter in the interval {a, a + da). Similarly, 
given a quantity aj that depends on the species index one replaces 



We next consider the possible methods for predicting the thermodynamic behavior of 
SHS fiuids. As pointed out in the introduction, a good approximation to the effectively 
exact Monte Carlo (MC) equation of state^° of the monodisperse SHS model is obtained 
by calculating the pressure from the energy route within the PY approximation^. In the 
case of mixtures no comparable Monte Carlo data exists, nor is a direct solution of the PY 
closure feasible, so that finding a reliable approximation to the equation of state remains an 
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important open challenge. As described in the introduction, we have tackled this in past work 
within an approximate theory based on a density expansion of the direct correlation function 
around the MSA solutio n^^i^^i^^ . Another possible route is thermodynamic perturbation 
theory. For the Baxter SHS model it is easy to convince oneself that only the scheme 
proposed by Weeks, Chandler and Anderson (WCA)^- can be applied. We have explored 
this possibility in the monodisperse case, where Monte Carlo simulations provide reliable 
reference data. In Fig. Q we compare the simulation data with the predictions of the Mean- 
Spherical- Approximation (MSA), the modified Mean- Spherical- Approximation (mMSA) and 
the CI approximation (as discussed inM-); the results from the first and second order WCA^^ 
perturbation theory are also shown. It is clear that the mMSA and CI approximations are 
fairly reliable for low and intermediate densities, even at low reduced temperatures, while 
the second-order WCA approximation breaks down already at temperatures significantly 
above the critical point (tc ~ 0.11, depending on the approximation used). The WCA 
method therefore offers little hope of providing the basis for an accurate equation of state 
for mixtures. One also sees readily from Fig. ^ that the PY closure provides by far the most 
accurate of all the approximation methods. This is why we return to the problem of solving 
the PY approximation for SHS mixtures in this paper. 

A major challenge in calculating phase equilibria in polydisperse SHS, or indeed any 
polydisperse fluid, arises from the fact that its Helmholtz free energy is a functional of 
the distribution p{a) of the polydisperse attribute^. However, in simple systems or ap- 
proximations this functional dependence reduces, for the excess free energy, to one on a 
finite number of moments of the distribution. In these cases the free energy is called trun- 
catable^Si^i and the phase coexistence problem reduces to the solution of a finite number of 
coupled nonlinear equations. For example, for the size-polydisperse SHS mixture the mMSA 
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and CI approximations yield such a truncatable form for the excess free energy involving 
only three moments pi,p2 and pa, and the two-phase coexistence problem can easily be 
solved numerically^^. The relevant moments are defined here including factors of density as 



Jo 

for m = 1, 2, 3; for later reference we note that ps is proportional to the hard sphere volume 
fraction. 

When the more accurate PY approximation is used, the presence of polydispersity renders 
an analytical calculation of the free energy impossible (see section HVBjl . In addition, even 
if the free energy could be calculated in closed form, it would almost certainly not have a 
truncatable form and so predictions for the phase behavior would remain difficult to extract. 
We therefore propose to consider a small degree of polydispersity as a perturbation^^ around 
the well-understood monodisperse reference system (see^^ for an overview of earlier work in 
this perturbative spirit). Denote by ctq a characteristic sphere diameter, which will be taken 
as the mean diameter of the overall or "parent" size distribution p^^\a) in the system. We 
then focus on fluids with a narrow size distribution centered on ao, for which the relative 
particle size deviations 



are small for all particle sizes a. Following Evans, we will expand up to second order in these 
size deviationsiS. The leading order phase boundary shifts and fractionation effects then turn 
out to be proportional to s^, where s = [(5^)'^°^]^/^ is the normalized standard deviation - 
also referred to simply as "polydispersity" - of the parent distribution. Before proceeding 
to the calculation, we address in the next section the choice of the stickiness coefficients e^j 
from Eq. These are irrelevant for monodisperse SHS but can have important effects on 




(3) 
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the behavior of mixtures as we will see. 

III. THE STICKINESS COEFFICIENTS 
A. General arguments 

At a reduced temperature r the Boltzmann factor exp[— 0jj(r)] for the interaction of 
two SHS particles depends only on the ratio ej^/r (and, of course, on cTjj). Physically, the 
stickiness coefficients represent dimensionless adhesion energies between pairs of particles 
identified by the species indices i and j. (We revert to the notation for the discrete mixture 
here; the same considerations obviously apply to the poly disperse system.) The e^j have no 
analogue in the monodisperse case, where only the reduced temperature r features and e 
can be set to unity. For (discrete or polydisperse) mixtures, on the other hand, one needs to 
make an appropriate choice for the dependence = jF(aj,crj) of the stickiness coefficients 
on the particle sizes. We discuss possibilities for this choice in this section. 

Clearly the appropriate form of the function T[pi^ aj) will depend on the kind of physical 
problem one is studying. Nevertheless, it should satisfy some general requirements: (i) 
Adhesion should be a purely pairwise property, and so JF should depend only on cij and 
aj as anticipated by our notation; JF must clearly also be symmetric under interchange of 
(Tj and (Tj. (ii) Since the eij are dimensionless, so must JF be. If it does not contain a 
separate lengthscale, it is therefore a homogeneous function of degree zero in {ai,(7j). The 
latter case is interesting because it can be seen as the sticky limit of a scalable (i.e. purely 
size-polydisperse) interaction^^, where by definition 4>ij{r) remains unchanged when r, cTj 
and aj are all scaled by a common factor. (The square well potential of Eq. can be put 
into this form by choosing Rij = <Jij[l + l/{Aeij — 1)]; the sticky limit is obtained by letting 
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A ^ oo.) The presence of pure size-polydispersity has important simphfying effects on the 
phase behavio r^^i^^ which we discuss further in Sec. below, (iii) If the adhesion depends 
on the surface area of the spheres one might expect JF to depend on ratios of homogeneous 
functions of degree two in (cTj, aj). (iv) If the adhesive interaction vanishes when at least one 
of the two particles i and j degenerates to a point we need to require limo-.^o -^(cTi, CTj) = 0; 
the limit for aj ^ is then also zero, by the symmetry of JF. 

In Ref.— plausibility and convenience arguments were adduced to suggest the following 
choices for the quantities eij: 

ol/afj Case I , 

Giaj/afj Case II , 
J^{au o-j) = { (5) 

1 Case IV , 

o^joij Case V . 

Here ctq is a characteristic reference length for the sizes, taken as above to be the parental 
mean diameter. In the forms originally suggested^^, this length was chosen as a moment of 
the size distribution, (cr")^/"' with either = 1 or 2. (Case I here corresponds to cases I and 
III in Ref."i^; we have kept the original numbering for the remaining cases II, IV and V for 
ease of reference.) However, this identification has the drawback of introducing many-body 
effects into the pair potential, as the moments (a") depend upon the thermodynamic state of 
the fluid, and in particular on the concentrations of all particle species. This is why we have 
chosen the fixed reference length above, consistent with the notion of a purely pairwise 
interaction. Numerically the actual choice of ctq turns out to have only very minor effects; 
this can be shown by calculations (not reproduced here) comparing case I (with fixed ctq) 
with case III from^^, obtained by replacing ctq — > (a^)^/^. 

The form of the denominator for cases I and II in Eq. ^ is forced by technical 
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constraints detailed in Ref.— , but these still leave some flexibility in the choice of numerator; 
cases I and II assume respectively a mean-fleld-like and a decoupled dependence between 
stickiness and size. Case IV corresponds to the choice of constant coefficients (independent 
of particle sizes), while case V was selected in Ref.'J^ specifically to permit analytical solution 
within the CI closure. Note that not all four cases have all of the properties (ii-iv) listed 
above as possible requirements. E.g. only cases II and IV are homogeneous functions of 
((Tj,(Tj) as required by (ii) when no additional lengthscale such as cxo is involved; they are 
therefore purely size-polydisperse. The properties (iii) and (iv) hold only for case II. It can 
be argued^i that the dependence on aiaj/afj assumed in case II is quite generic for solutions 
of colloids, micelles or globular proteins, at least in the high-temperature regime where a 
linearized approximation for the Boltzmann factor is sufficient. While this favors case II, for 
phase coexistence we are interested in lower temperatures where it is less clear which case 
is physically more appropriate; we will therefore include all four cases in our analysis. 

For our perturbative analysis we only need to know the coefficients in the expansion of 
the eij around the typical particle size o"j = aj = ctq, up to quadratic order in the relative 
particle size deviations 5i = (aj — ao)/aQ: 

(^ij = Co + eia(5i + Sj) + e2a5i5j + e2bi5i + 5|) + . . . (6) 

The coefficients eo, eia, e2a, ^2b of this expansion are given in Table IHfor the four cases listed 
above. Note that eo = 1 always so that in the monodisperse limit the e^j are irrelevant as 
they should be. 
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B. Stickiness coefficients for the Asakura-Oosawa model 



So far we have considered choices for the stickiness coefficients suggested by rather general 
arguments. One may wonder whether the e^j can be derived more directly from a physical 
picture. We shall pursue this here for the well-known Asakura-Oosawa model of colloid- 
polymer mixtures, which for small polymers leads to a short-ranged attractive depletion 
potential acting between the colloids^. We shall show that, while a formal sticky limit 
cannot be taken in general when colloids of different sizes are present, an effective SHS 
model can still be derived when the polymer size is small but kept nonzero. This then 
simplifies further in the perturbative approach for weak polydispersity adopted here. 

Consider two colloidal particles represented by impenetrable spheres of diameter cTj and 
aj immersed in a solution of non-interacting polymers. Within the Asakura-Oosawa model, 
the polymers are simplified to spheres of diameter ^ which can fully penetrate each other but 
have a hard sphere interaction with the colloids. It is well known that such a system develops 
an entropically driven effective attraction between the colloidal particles. This arises due to 
a reduction in the volume from which the polymers are excluded when the exclusion zones 
around the colloids overlap (see Fig. |2)). This overlap volume as a function of the distance 
r between the sphere centers is 



where Rk = {ck + 0/2 and only distances r > 0"^ are allowed because of the hard colloid- 
colloid repulsion. The effective colloid-colloid attraction induced by the presence of the 
polymers is then just the overlap volume times the polymer osmotic pressure^^^, giving the 



Vov(r) 



12 



TT 



- ^{R\ + R^^r + %{Pi\ + - 3(i?,' - R^^f- e{aij + ^ - r) (7) 
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(8) 



overall AO interaction potential 

+00 < r < (Jij , 

-PpVov(r) (Jjj < r < (7ij- + ^ , 
r > (7jj + ^ , 

This expression can be obtained formally by integrating out the polymer degrees of freedom 
from the partition function at fixed polymer chemical potential. The latter is conveniently 
parametrized by the density Pp of polymers in a reservoir connected to the system; because 
the polymers are taken as ideal, their osmotic pressure is then k-oTpp and the U-qT is absorbed 
by our choice of units. The effective colloid-colloid interaction will in general contain also 
many-body terms, but these vanish in the limit of small polymers (for monodisperse colloids 
the condition is ^ < 0.1547(Jo) that we are interested in. 

To map to an equivalent SHS potential, which should be physically reasonable for small 
polymer-to-colloid size ratio ^/ctq, one equates the corresponding second virial coefficients. 
The hard core makes the same contribution {B2,'rs — '^'^^Ij/'^) SHS and the original 

AO potential, so one can focus on the normalized deviation of the second virial coefficient 
from this HS value. 



^-°2,AO 



-"2,AO -"2,HS 

'13 '' <^ij 



For the SHS potential this quantity equals — l/(4rjj), so the stickiness parameters in the 
mapped SHS system are assigned as 

1 1 





-4,^0 (r) ^ 


dr . 


(9) 


1 aij 









We now proceed to simplify this expression for small ^; in the limit ^ ^ 0, the original 
AO model should become fully equivalent to the mapped SHS system. We will see that 
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for mixtures of colloids of different sizes this strict mathematical limit cannot be taken 
consistently; nevertheless, as long as ^/ctq is small, we expect the SHS mixture to give a 
reasonably accurate description of the underlying AO model. 

To simplify Eq. we change integration variable from r to z = {r ~ aij)/^, expand the 
attractive tail of the AO potential in ^ as 



ij 



a, 



z-if+o{e) ■ 



(10) 



and retain only the leading term. Similarly approximating = [aij + ^z)'^ = afj + 0{^) 



yields 



12r,;, 





f 




dz = — 


CTij 
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where 



is the value of the attractive potential at contact and erfi(z) = erf {iz)/i is the imaginary 

error function. Because of the prefactor ^/(Jij in Eq. (fTT|) . 'jij has to grow as C, decreases 

if we want to keep r^j finite. For large argument the error function behaves as erfi(z) = 

e'^[l/z + 0(1/ z^)]/^ and so 

1 e e^- 2 
= e »j . 

A nonzero limit value of Tij for ^ ^ thus requires that jij grows logarithmically as 'jij = 
ln{aij/C,) to leading order. The corresponding polymer reservoir density, likewise to leading 
order, goes as 

4 aij ln(cro/0 



Pp 



(12) 



The dominant dependence Pp oc C, ^ in this expression arises because the value of the AO 
potential at contact scales as Pp^"^] the additional logarithmic factor increases this interaction 
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strength to compensate for the decreasing range of the attraction as ^ 0. Note that even 
though the polymer density diverges, the polymers do in fact become very dilute as one sees 
from the (reservoir) volume fraction rjp = {'n'/Q)pp^'^ ~ .^ln(cro/.^) occupied by the polymer 
spheres. 

For monodisperse colloids, the above procedure produces an unambiguous sticky limit 
for ^ ^ 0. The explicit form of Eq. (fT^ shows, however, that this limit cannot be taken 
straightforwardly for mixtures: the prefactors cTjj/ {o-icrj) of the required leading order diver- 
gences of the polymer density are incompatible with each other for different pairs of particle 
species. In other words, if the ^-dependence of the polymer density is chosen to keep one 
specific Tij finite and nonzero, then the others would either tend to zero or grow to infinity 
in the sticky limit. The example of a binary mixture illustrates this. Suppose that ai > 02 
and that the polymer density is tuned to keep the rn finite. Then l/ri2 and l/r22 would 
both tend to zero for ,^ ^ so that all interactions involving particles of species 2 become 
purely HS-like, without any attractive contributions (this is system B studied in2^). 

In the absence of a strict sticky limit, we will content ourselves with applying the mapping 
|TT|) for small but nonzero polymer-to-colloid size ratios ^/ctq. The properties of the resulting 
SHS mixture should then still give a good approximation to those of the original AO model. 
In the perturbative setting of this paper we can then expand Eq. (fTT|) in the small relative 
deviations 8i = (crj— (Tq) / o"o of the particle sizes from the parental mean. In the decomposition 
1/Tij = etj/r of Eq. (jSJ we fix the scale of the e^j by requiring as before that e^j = 1 for 
particles of the reference size = aj = ctq. This gives 




(13) 
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for the reduced temperature, where 

7 = -^Pp^ (^0 ■ 

The second, approximate equahty in Eq. ()13|) holds for large 7 as before. To find the 
perturbative expansion of the stickiness coefficients tij, we note first that the potentials at 
contact expand as 



1^J = 7 



2' 2 ' 4 

Since the erfi in Eq. (jlll) grows at most as exp(7jj), a second order Taylor expansion will 

give an accurate approximation as long as the perturbations in jij are -C 1. This requires 

5j ^ 1 /7, which then automatically enforces 6i <^ 1 since we expect 7 to be at least of order 

unity for the mapping to a SHS mixture to make sense. Under these conditions one then 

has a valid perturbation expansion of the e^j. The coefficients defined in Eq. (0) are found 

as eo = 1 (by our choice of r) and 

_ -1+ gi _l + 92 _l-2gi+ g2 
Ci — ; ^2a — — — ? ^2b — : 



-1 1 

91 

TT 



where 



/7erfi(^)-2 2 
_ [3 + e^(27-3)]/4 3 

92 — — r ~ . 

A/7r/7erfi(y7) - 2 8 
From Eq. (fT^ one sees that the reduced temperature is set by the contact potential 7, which 
itself is proportional to the polymer reservoir density. Unlike the more ad-hoc choices of 
Eq. (0), the expansion of the in terms of the bi depends on the reduced temperature r, via 
7. For large 7 one can use the leading order approximations (71 ^ 7 — 1, (72 ~ (7^ ~ 27 + 1) /2 
to evaluate this dependence. However, since typical values of 7 are only logarithmically 
large in o^j^ it is generally safer to work with the full expressions. 
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IV. PERTURBATION THEORY FOR THE POLYDISPERSE PY CLOSURE 

In this section we come to the core of our analysis. We first review Evans' perturbative 
framework for shghtly polydisperse systems. To apply this to the PY approximation for 
SHS mixtures we will need the perturbative expansion of certain correlation function values 
at contact; from these we can then finally find the excess free energy. 

A. Evans' perturbative expansion 

The starting point for an analysis of the phase behavior of polydisperse systems is the 
excess free energy density. In general this is a functional of the size distribution p{a) in the 
system. It is also a function of the particle density p, and of temperature; we do not write 
the latter explicitly below. For slightly polydisperse systems it is expedient to switch from 
a to the relative deviations 6 from the reference size (Tq. By the fundamental assumption 
of a narrow size distribution, the 5 are small quantities, and one can expand the excess free 
energy density f^^, measured again in units of k^T, in terms of moments of 

rip, [pm = rip) + paip){6) + pbip){6') + pcip){5f + . . . . (14) 

Here terms up to second order in 5 have been retained; these give the leading effects on 
the phase boundaries'^. Our functions a, 6, c differ by factors of p from those defined in 
Ref.— , so that e.g. a equals Evans' A/ p] this simplifies the statement of Eqs. ()15II17|) below. 
The leading term /g^ is the excess free energy density of the monodisperse reference system 
where all particles have 5 = 0. 

Given the above expansion of the excess free energy, the conditions for two-phase equilib- 
ria of the near-monodisperse fluid can be solved perturbatively^S,. We briefly recall the main 
results. The fluid is initially in a parent phase of density p*^°^ , with a parent size distribution 
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function p*^") (5) , where (5) = by our choice of the reference size ctq as the parental mean. 



In order to lower its free energy, the fluid can split into two daughter phases of densities p' 



,(1) 



and p*-^-* , with distribution functions p*-^-* (5) and p*-^-* (5) which are in general different from 
the parent distribution, a phenomenon referred to as fractionation^^ . The densities and size 
distributions can be worked out perturbatively at any point inside the coexistence region^^; 
we focus on the properties at the onset of phase coexistence, which are most easily accessible 
experiment ally. 

Suppose the system is just starting to phase separate, with all of the volume except for an 
infinitesimal fraction still occupied by phase 1, with density p^^\ Conservation of particle 
number then requires that p^^\5) = p^^\6), i.e. the size distribution in this cloud phase 
equals the parent. The coexisting shadow phase 2, on the other hand, will generally have 
7^ p^^\6). Evans2^ showed that the cloud and shadow densities, p'^^^ = p^^^ + 6p''^^ 
are shifted from their monodisperse values p^^^ and p^^^ by 



Here a' = da/dp, b' = db/dp and /t(p) = l/[p + p'^{d/dpy f^^^p)] is the isothermal compress- 
ibility of the monodisperse reference system. The shorthand A indicates differences between 
the two monodisperse reference phases, e.g. Aa = a(po^^) — a{p^Q^). Finally, recall that s is 
the parent polydispersity: the phase boundary shifts are to leading order quadratic in s. 

It is worth noting that Eqs. are not symmetric in p^Q^ and Po^^; by interchang- 

ing the two densities one therefore obtains a different cloud-shadow pair. Physically, this 
corresponds to approaching the onset of phase separation from low or high densities; in a 
polydisperse system the coexisting phases are different in the two situations since only the 
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respective majority (cloud) phase has the parental size distribution. The size distribution 
in the corresponding shadow reads, to leading order in 6^, 

=p(o)(5)[i + (Aa)(5] . (17) 

Overall, the monodisperse binodal delimiting the coexistence region splits into separate 
cloud and shadow curves, which intersect in the critical point Quantitative information 
about the critical region is not accessible within the perturbative expansion of Eqs. ()15ll6p . 
however, since the compressibility k diverges as the critical point is approached. 

The above summary shows that knowledge of the functions a, b and c is sufficient to 
calculate the leading order phase boundary shifts and fractionation effects for weakly poly- 
disperse systems. In the next two subsections we calculate these functions for the SHS 
mixture within the PY approximation. 

B. Perturbative analysis of the PY closure 

To lighten the notation in the rest of the paper, we make all densities dimensionless by 
measuring them in units of Vq^, where 

is the volume of a particle with the reference diameter. The third moment ps defined in 
Eq. is then identical to the hard sphere volume fraction rj. We also measure all particle 
sizes a in terms of ctq, so that the relation between a and the fractional deviation from 
the parental mean diameter becomes simply a = 1 + 6. In the monodisperse case, where 
all particles have 6 = 0, all moments (jH)) are then identical and equal to the density p 
(which also equals the volume fraction rj). Finally, for notational simplicity we again revert 
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temporarily to the case of a discrete p-component SHS mixture; the final results will be 
expressed in terms of averages over the size distribution and so generalize immediately to 
fully polydisperse systems. 

In order to extract the desired thermodynamic quantities from the PY closure, the fol- 
lowing set of p{p + l)/2 coupled quadratic equations needs to be solved firsli^, 

Lij ttjj -|- ^ ^ . 

m 

where the unknowns are 



Xri 



m 



z,j = l,2,...,p (18) 



Lij 



T 



Here Vijio'ij) is the partial cavity function at contact which is proportional to the probability 
of finding a particle of species j touching any given particle of species i. In Eq. (jl8|) the 
coefficients aij, jSij, and (pij are given by 

Pij = puijeij/r , (20) 
4>ij = <^i<^j/^ ■ (21) 



Here the quantities 



are the PY partial cavity functions at contact for the HS fiuid (to which the SHS fiuid 
reduces at infinite reduced temperature r) and we abbreviate A = 1 — 77, with r] = the 
HS packing fraction as before. Notice that all four sets of coefficients Lij,aij, f3ij, and (pij 
are symmetric under exchange of the species indices i and j. 

For one-component fiuids, the system (fTHj) reduces to a single quadratic equation. Baxtei^ 
showed that only the smaller of the two real solutions (provided such solutions exist at all) 
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is physically significant; it is given explicitly in Eq. ()24p below. For true mixtures {p > 1), 
explicit solution of the rather complicated system (fTHj) of algebraic equations is feasible at 
best numerically (except for special casea^^^) and is the computational bottleneck of the 
PY solution. For large p, and certainly for the polydisperse limit p — > oo, it is impossible 
in practice. However, progress can be made for near-monodisperse fluids by solving (fTH)) 
perturbatively. The Lij will generically depend on the reduced temperature r, the overall 
number density p, the sizes cTj and aj of the particles at contact, and all the molar fractions 
Xi (or their polydisperse analogue, the size distribution p{S)). For small Si we can therefore 
expand to quadratic order as 



The idea now is to insert this expansion, and the analogous expansions of the known co- 
efficients aij, Pij and (f)ij, into the r.h.s. of Eq. (fTSj) . Having done this, one re-expands 
to quadratic order in 6i, 6j, 6m and (6), and to linear order in (5^). Finally one replaces 
^m^rn = 1 and 'Ylm^m^m = (f^") fo^ = 1)2. Comparing terms of the same form on 
the left and right of Eq. (fTH|) one then finds a relatively simple set of equations for the 
coefficients Lq, . . . , L2e, as outhned in the Appendix. To order zero in polydispersity one of 
course retrieves Baxter's original quadratic equation (Eq. ()A5|) ). whose physically relevant 
solution is 



where Aq = 1 — p is the value of A in a monodisperse system with density p. Since we are 
perturbing around the physical solution for the monodisperse case, the results we find 
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for slightly polydisperse mixtures will automatically have the correct physical behavior. In 
a non-perturbative solution, one would need to check separately that the solution branch 
with the correct low-density limit Lij —>■ ufj/rij has been selected; this condition arises since 
Viji^ij) ^ 1 at low density. 

The conditions imposed by Eq. ()18|) for the higher order expansion coefficients Lia, . . . , 
turn out to be linear and can be straightforwardly solved order by order; see the Ap- 
pendix. The region in the density-temperature plane where Eq. (|TH|l has no physical solution 
therefore remains as in the monodisperse case, being delimited by < p < with 



= 5-12r + 6r^ " ^'^^ 

for r < (2 — v^)/6. This is clearly an artifact of our finite-order perturbation theory, given 
that we know from numerical solutions of Eq. ()18|) that the region where solutions exist 
does change with increasing polydispersityi^. To reproduce this effect within our approach, 
a resummation of the perturbation theory to all orders would be needed. 

C. Excess free energy 

Given the perturbative expansion for Lij, we can determine the free energy of weakly 
polydisperse SHS mixtures in the PY approximation. There are three known thermodynamic 
routes (via the energy, compressibility, and virial) that could potentially be used^i. We focus 
on the one that gives the most reliable equation of state for the monodisperse system (see 
Fig. IH), i.e. the energy route. It predicts in general for the r-derivative of the excess free 
energy density 
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Inserting the expansion (j23jl of Lij and re-expanding to quadratic order yields 

or 



dr 



where 



To — Lq 



Ti = Lo + 2Li„ + L 



lb 



^2 — Lia + Lib + + 2L2c + -^^2(i , 



Ts — + 2L26 + 



We can then integrate from the desired value of r to the hard-sphere limit r — 00 to find 



where 



A/r = V/ r.(rO^ 



i = 0,1,2,3 . 



and /^g is the excess free energy density of the HS fluid. For the latter we use the standard 
Boublik^, Mansoori, Carnahan, Starling, and Leland!2^ (BMCSL) expression^. Expanded 
to second order in polydispersity this reads 

/hs = /hs,o + /hs,i W + fnsM + /hs,3('^') , 
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Altogether we therefore have, for the perturbative expansion ()14|) of the excess free energy 
density of the SHS mixture, 

/ex xex I A xex 

Jo — iHS,0 "T ^iO ; 

ap = f^s,i + A/r , , , 

(26) 

bp = /fg,3 + A/r , 

cp = f^s^, + A/r . 

With these results we can now proceed to apply Evans' general results to study cloud and 
shadow curves and fractionation effects in polydisperse SHS mixtures. 

Inspection of the lengthy explicit expressions for a, b and c shows that the dependence 
on the stickiness expansion coefficients eia, e2a, ^2b is in fact rather simple. For a one finds 
the form 

a = ao + eiatti , (27) 

with qq and ai functions of p and r only. This is reasonable since a is the coefficient of a 
first order (in S) term in the excess free energy, and should therefore only depend on the 
expansion of the tij to the same order. The function b involves in addition terms proportional 
to el^ and e2b, while the remaining coefficient occurs only in the function c. Since c does 
not feature in the expressions for the phase boundary shifts or fractionation effects to O(s^), 
all results we show below are therefore independent of e2a- 



V. PHASE BEHAVIOR 



In this section we show our results for the phase behavior of polydisperse SHS mixtures. 
We will explore the various choices of stickiness coefficients discussed in Sec. IIIH i.e. cases 
I, II, IV and V as well as the AO model for small values of the polymer-to-colloid size ratio. 
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The first subsection has the main results from our perturbation theory in polydispersity 
for the PY closure; in Sec. IV Bl we then compare these predictions with those from other 
approximation schemes. 

A. PY closure 

We start by recalling in Fig. IHl the phase diagram of the monodisperse SHS fluid as 
obtained within the PY approximation and using the energy route to thermodynamics. 
Along with the binodal we show the spinodal, where the curvature of the free energy vanishes 
and a homogeneous phase becomes unstable to local density fluctuations, and the region p5|l 
where Baxter's PY equation has no physical solution. Here and in the following we use 
on the X-axis the volume fraction rj rather than the density p. In our units, these two 
quantities are identical for monodisperse systems, but differ to order in the presence of 
size polydispersity. For parent phases specifically, Eq. ()A4|) gives 77*^°^ = p*^°^(l + 3s^) to 
quadratic order. Cloud phases, which share the parental size distribution, have similarly 
P^^^ = Po^\l + 3s^) + Sp^^\ while for shadow phases one finds using Eq. (|T7jl that p*^^^ = 

p!,'^[l + 3(l + Aa)s2]+5p('^^- 

To get some initial intuition for the effects of polydispersity, it is useful to consider 
first the single-phase equation of state. Fig. |^ shows plots of the dimensionless pressure 
against volume fraction at several values of the polydispersity and for three choices of the 
reduced temperature r. We consider here constant stickiness coefficients (case IV) to allow 
a comparison with numerical work for discrete mixtures^^,. It is gratifying that we find 
quahtatively the same trend, with the pressure decreasing with increasing polydispersity. 
Quantitatively, however, the results are not directly comparable because in Ref.~ the less 
accurate compressibility (rather than energy) route was used to evaluate the pressure. 



26 

To interpret physically why the pressure decreases with polydispersity s at fixed packing 
fraction rj, we note first that such a decrease is found also in the absence of adhesion (i.e. 
for HS). This has been established in simulationa^ and is reproduced qualitatively by the 
BMCSL equation of state; the intuitive reason is that in a fiuid (gas or liquid) phase a 
spread of sizes allows for a more efficient packing of the particles. In such a less "jammed" 
particle arrangement one expects to find fewer interparticle contacts and so, in the presence 
of adhesion, fewer particle pairs interacting attractively. This will increase the pressure, 
counteracting the reduction, that one would expect for HS, resulting from the more efficient 
packing. Our results are quite consistent with this: at finite r, we find that the pressure 
decreases less with polydispersity than in the HS limit r — * oo. 

The curves shown for the polydisperse cases in Fig. 0] cannot be used to infer phase 
coexistence properties directly by e.g. a Maxwell-construction: fractionation means that two 
coexisting phases do not have properties represented by a single relation between pressure 
and volume fraction. This remark holds true quite generally for single-phase equations of 
state in polydisperse systems, including e.g. the results obtained in Ref.^^ within the PY 
compressibility route to the equation of state. However, some more limited information on 
single-phase stability can be deduced. Specifically, a single phase cannot be stable where the 
pressure decreases with volume fraction. For the middle graph of Fig. |5 for example, where 
T = 0.1186 is just above the monodisperse critical point and so a monodisperse system is still 
stable at all densities, the polydisperse mixtures with s = 0.2 and 0.3 are already unstable 
in some range of densities. This means that the region where phase separation occurs must 
extend to larger values of r for polydisperse than for monodisperse SHS, a result which - 
for case IV, as considered here - we will find confirmed very shortly. 

We next turn to explicit results for the phase behaviour, starting in Fig. El with cases 
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II and IV for the stickiness coefficients, illustrated here for parent polydispersity s = 0.3. 
The cloud curve gives the boundary of the region where phase coexistence occurs. The 
shadow curve, which records the density of the coexisting phase at each point of incipient 
phase separation, is normally distinct from this. However, for the purely size-polydisperse 
cases considered here it is known on general grounds that when represented in terms of 
volume fraction rather than density the cloud and shadow curves coincide to O(s^)'^^'^'^. It 
is reassuring that, as Fig. shows, this property is preserved by the PY approximation. 

Turning to more detailed features of Fig. El we observe that in case IV the coexistence 
region is broadened towards both lower and higher volume fractions. As the monodisperse 
critical point is approached, the perturbation expansion breaks down as expected and the 
cloud/shadow curves diverge. No quantitative information can then be extracted in this 
regime, but the fact that the divergence is outwards still tells us that the coexistence region 
in the polydisperse case extends to larger values of r than for monodisperse SHS. This is 
consistent with our inference from the single-phase equation of state above. 

Comparing cases II and IV in Fig. El one sees first that the phase boundary shifts are 
rather smaller in the former than the latter. Also the (slight) broadening of the phase 
separation region towards lower rj is now restricted to r below around 0.093, while above 
the opposite trend is observed. The divergence of the curves at the monodisperse critical 
point is now inwards so that phase coexistence must terminate at a values of r below the 
monodisperse Tc- 

Figure El shows the cloud and shadow curves for case V. We find that the shifts away from 
the monodisperse binodal are rather larger than in the previous two cases, and therefore 
show results for a smaller polydispersity s = 0.2, rather than for s = 0.3. Cloud and 
shadow curves no longer collapse, consistent with expectation as case V is not purely size- 
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polydisperse. The cloud curve shows that the coexistence region narrows in this case, except 
on the high-density branch for r below ^ 0.085. The inward divergence of the cloud curve 
shows that the coexistence region also shrinks towards lower r. The shadow phases are more 
dense throughout than the phases on the same branch of the cloud curve. Except for the 
last point, these trends agree with the non-perturbative results of Ref.— derived within the 
CO closure. 

Case I, shown in Fig. [3 has even stronger polydispersity effects and we show predictions 
for a correspondingly smaller polydispersity s = 0.1. For r not too far below the critical 
point the behaviour is otherwise qualitatively similar to case V; for lower r the coexistence 
region is displaced towards lower rather than, as in case V, higher volume fractions. The 
shrinking of the coexistence region towards lower r is again in qualitative agreement with 
results from the simpler CO closure^^. 

Finally we turn to the phase behaviour predicted for the AO model with a small polymer- 
to-colloid size ratio ^/o"o = 0-1 ^'^'^ polydispersity s = 0.07, as shown in Fig. |H| For this 
choice of ^ we have 7 ~ 3.97 at the critical point of the monodisperse system, and the 
condition (5j ~ s ^ I/7 for the validity of the expansion in particle size of the stickiness 
coefficients e^j is reasonably well obeyed. Here the coexistence region is broadened in all 
directions by the introduction of polydispersity: towards low and high densities, and also 
towards larger values of r. The shadow phases are again more densely packed than the 
analogous cloud phases. 

We conclude this section by considering fractionation effects. These are illustrated in 
Fig. El for cases II and I, for a parent distribution of Schulz form and with values of the 
polydispersity s as in the corresponding Figs.Eland[7| When phase separation is approached 
from low densities, a gas cloud phase with the parental size distribution coexists with an 
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infinitesimal amount of a liquid shadow phase with a different size distribution. At the 
high density boundary of the coexistence region, a liquid cloud phase similarly coexists with 
a distinct gas shadow phase. Fig. IHl shows that for case II the liquid phase contains more 
larger particles than the coexisting gas in both of these situations (and therefore presumably 
throughout the whole coexistence range of parent densities at the chosen r). Case I exhibits 
the opposite behaviour: here the liquid phases contain more smaller particles than their 
coexisting gas counterparts. 

To understand this difference between cases I and II, we return to Eq. (fT7j) . Consider 
the gas cloud point, where Pq^"* and p^^^ are the densities of coexisting gas and liquid in the 
monodisperse system; Aa then is the difference in the values of a between gas and liquid. 
If this is positive, then Eq. (|T7jl says that the liquid shadow has an enhanced concentration 
of larger particles. By reversing the role of the two densities one then sees easily that also 
at the liquid cloud point the liquid phase will contain more of the larger particles than the 
gas (shadow) phase. In summary, the liquid contains predominantly the larger particles if 
Aa > 0, and the smaller particles if Aa < 0. But from Eq. (jTTj) . Aa = Aao + eiaAai so that 
different choices of stickiness coefficients affect the direction and strength of fractionation 
only via ei^. The functions Aag and Aai are shown in Fig. ^Jand are both positive; as a 
result, Aa is positive when eia > — Aao/Aai and negative otherwise. The ratio occurring 
on the r.h.s. is almost constant and remains close to 1/3 over a large range of r, as the 
inset of Fig. El demonstrates. We can now rationalize the difference between cases I and II 
observed above: for case I, eia = — 1 < —1/3, hence Aa < and fractionation will enrich 
the liquid in small particles; for case II, ei^ = > —1/3 and one has the opposite situation. 
Referring to Table H] we also conclude that case IV will have the same fractionation behaviour 
as case II, while case V will produce the same "direction" of fractionation (smaller particles 
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in the liquid) I but with quantitatively weaker effects. In the AO case eia depends 

on r as discussed in Sec. IIIIBl but this effect turns out to be weak quantitatively, with (for 
^/(jQ = 0.1) eia ranging from ^ 0.95 at the critical point to ~ 1.24 at r = 0.065. Taking for 
simplicity eia ~ 1 one infers that fractionation effects will be qualitatively similar to cases II 
and IV, but quantitatively Aa will be larger by a factor of around 4. All of these conclusions 
can be confirmed by detailed examination of the explicit results for the various cases. 

B. Other approximation schemes 

Once one accepts the PY closure, the results shown above are exact in their treatment 
of polydispersity, certainly within the perturbative setting of weakly polydisperse mixtures. 
However, the PY closure itself - while more accurate than its competitors - does remains an 
approximation. It is therefore useful to compare with the predictions of other approximation 
schemes to assess the robustness of our predictions. We do this first for case II, where an 
approximate free energy of BMCSL type can be constructed, and then for the AO model, 
which can be analysed using the free volume theory of Refs.-^^^. 

To construct the alternative approximation for case II one starts from a virial expansion 
of the excess free energy density up to the third virial coefficient. This is easily found as 

r = pp3 + (3 - i2t)pip2 + 

^ [ppI + 3(1 - 12t + 48^2 - 32t')pl + 6(1 - 4t)pip2P3] , (28) 

where t = 1/ (12r); the terms of second order in density agree with the energy route of the 
CO approximationi^. The interesting feature of this result is that the fourth order moment 
P4 does not appear, in contrast to the analogous expansions for the other cases I, IV and 
V that we have considered. Furthermore, the only modification compared to the pure HS 
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case is in the t-dependence of the coefficients. These observations suggest that it should 
be possible to construct a modified free energy expression of BMCSL-type which matches 
the above virial expansion to third order in density. Remarkably, if the desired modified 
BMCSL form is parametrized in a fairly general manner as 



r ^ ( - M,) „„,1 - O,,) ,E].^^. ^^J^ . ,29, 



then by expanding to third order in density and matching to the expansion ()28p one finds a 
unique solution for the coefficients: 

^ = 0, D = A2 = l , B = l-At, C = Ai = + 32t^ . 

The presence of polydispersity is crucial here: for a monodisperse system, the matching 
conditions to third order in density would not constrain the coefficients sufficiently. 

One can now apply the perturbative scheme used throughout this paper to obtain from 
the excess free energy of Eq. fOHjl the functions a and b, and hence the cloud and shadow 
curves. (Note that the perturbative approach is used here mainly for ease of comparison with 
our other results; since the free energy ()29|) is truncatable, a full solution of the phase equi- 
librium conditions would be fairly straightforward.) The results are shown in Fig. ^2 note 
that not just the polydisperse cloud/shadow curves but also the monodispere binodal are 
different from the ones obtained from the PY approximation. Looking at the polydispersity- 
induced shifts, one sees that on the high-density branch of the cloud/shadow curve these are 
quite comparable to those from the PY approximation (Fig. EJ, even semi-quantitatively. 
Polydispersity effects on the low-density branch are rather smaller, again as found within 
the PY closure. Near the critical point, however, the trends are reversed: the BMCSL-type 
approximation predicts an extension of the coexistence region towards larger r and smaller 
Tj, whereas the PY approximation leads to the opposite result. 
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The second case where we have an alternative approximation scheme available for com- 
parison is the AO model. The free volume theory of Ref.~ effectively linearizes the excess 
free energy in the polymer (reservoir) potential pp, and the same is true for its generaliza- 
tion to polydisperse colloids^^. It is therefore most accurate when the depletion interaction 
between the colloids, which is proportional to Pp, is small (in units of k^T). In order to still 
get gas-liquid phase separation, the polymer size ^ must then not be too small. This is the 
opposite limit as for our SHS mapping, which will work best when ^ <^ ctq and the depletion 
attraction is large at contact. If anything one therefore expects the best agreement between 
the two approximations for intermediate values of ^; a suitable choice is ^/cxo = 0.1 as inves- 
tigated above. Fig. [T21 compares the two sets of cloud and shadow curves predicted. On the 
vertical axis we plot the polymer (reservoir) volume fraction rjp. This equals pp^^ in our di- 
mensionless units and is the conventional variable used in phase diagrams of colloid-polymer 
mixtures^^. Comparison of the two panels of Fig. 1121 reveals that the qualitative agreement 
between the two theories is surprisingly good. In particular, the qualitative changes caused 
by the presence of polydispersity (broadening of coexistence region to lower and higher col- 
loid volume fraction and lower polymer volume fraction) are in full agreement. For the 
relevant range of polymer volume fractions there is even quite good quantitative agreement 
(though note the slightly different axis ranges on left and right), and the shifts of cloud and 
shadow curves away from the monodisperse binodal are also quite comparable. Even the 
predicted fractionation effects agree well: as the inset on the right of Fig. [T21 demonstrates 
the calculated values of Aa are, apart from the slight shift in the critical point values of the 
polymer volume fraction, quite consistent with each other. 

We note briefly that in order to calculate the free volume theory data shown in Fig. 
we took the excess free energy for fully polydisperse colloids (at fixed polymer chemical 
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potential) derived in Ref.— and then found the functions a and h by expanding exphcitly 
as in Eq. ()14|) . This gives for a the same result as obtained by Evans^^, while h differs 
from his expression in terms of approximate correlation functional. One might expect that 
our approach of deriving a and h from one unified polydisperse excess free energy would be 
somewhat more accurate than Evans' procedure of finding a and h by quite different routes. 
We have checked that for larger polymer sizes ^/ctq = 0.4 our method predicts similar trends 
to those reported in^^, but quantitatively the effects of polydispersity are less pronounced. 

VI. CONCLUSIONS 

We have presented a perturbative approach to the determination of the gas-liquid phase 
behaviour of polydisperse Sticky Hard Spheres (SHS), studied within the Percus Yevick (PY) 
integral equation theory. For arbitrary size polydispersity, the calculation of phase diagrams 
analogous to those reported here would normally require the solution of a large (or infinite) 
system of quadratic coupled equations, a task which in practice can be accomplished neither 
analytically nor numerically. To get around this bottleneck of the PY closure we focussed 
on weakly polydisperse mixtures, where the overall size distribution is narrow in the sense 
that its normalized (by the mean) standard deviation s is small compared to unity. This 
allowed us to calculate in closed form the leading order (0(s^)) shifts of cloud and shadow 
curves away from the monodisperse binodal, and the corresponding fractionation effects. 
The thermodynamics was derived from the PY solution via the energy route because in the 
monodisperse case this method gives the best match to Monte Carlo simulation results, even 
for low reduced temperatures r around and below the critical point. 

In order to specify the properties of a SHS mixture one needs to know how the stickiness 
coefficients eij depend on the sizes of the two interacting particles. We discussed a number 
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of plausible constraints on this size dependence. In obtaining explicit results we considered 
specifically the cases I-V (excluding III, which with our now more appropriate choice of 
reference length becomes identical to I) previously suggested within exact solutions of simpler 
closures like CO and CI. Of these, cases II and IV are special since they can be seen as the 
sticky limit of purely size-polydisperse interactions, in which scaling of both particle sizes by 
a common factor only changes the range but not the strength of the interaction. We have 
also considered the AO model of a mixture of polydisperse colloids and polymers, which 
for small polymer size can be mapped to a good approximation onto an SHS model. The 
stickiness coefficients can be derived in this case rather than postulated; in contrast to the 
simpler ad hoc prescriptions of cases I-V, they are functions of r. 

In the simplest case IV of constant stickiness coefficients we first investigated the single- 
phase equation of state, finding qualitative agreement with a numerical solution of the 
compressibility equation of state for a small number of components by Robertus et al.— . 
Moving on to phase coexistence proper, we found for cases II and IV that cloud and shadow 
curves coincide in in the volume fraction representation and to O(s^), as expected on general 
grounds; less obviously our results also show that in these two cases the deviations of the 
polydisperse cloud/shadow curves away from the monodisperse binodal are quantitatively 
small. In all the other cases considered the shadow curves are located at higher volume 
fractions than the cloud curves, a trend observed in many other polydisperse system o^^i^^ . 

Summarizing our findings regarding the effect of polydispersity on the extent of the 
coexistence region as delimited by the cloud curve, it is simplest initially to group the 
different scenarios according to their behaviour near the critical point. For case IV and 
the AO model (with a polymer-to-colloid size ratio of 0.1) the coexistence region is shifted 
to higher reduced temperatures r; conversely, at fixed r it covers a wider range of parent 
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volume fractions 77. Cases I, II and V, on the other hand, show the opposite behaviour, with 
the coexistence region shrinking towards lower r. 

The trends in cases IV and AO remain unchanged as one moves to lower values of r, with 
the coexistence region continuing to broaden towards lower and higher values of rj at the 
two ends (gas and liquid). In the other cases the shrinking trend near the critical point can 
be reversed at lower r. E.g. for case II one also eventually sees a broadening to lower (gas 
branch) and higher (liquid branch) rj. For case V the coexistence region is shifted to higher 
rj at both ends (gas and liquid) at low r; case I shows the opposite behaviour. 

We have analyzed also the fractionation effects that accompany polydisperse phase sepa- 
ration, where coexisting phases have different particle size distributions. Depending on the 
stickiness coefficients considered, the liquid phase contains predominantly the larger (as in 
cases II, IV and AO) or the smaller particles (as in cases I and V). We rationalized this 
result by showing that the fractionation effects depend on the stickiness coefficients only via 
the expansion coefficient ei^; where this is above ~ —1/3, the larger particles accumulate in 
the liquid phase, otherwise in the gas. 

Finally we have compared our results with the predictions from other available approx- 
imation schemes, to check their robustness. Case II is important here because a variety 
of simple but realistic interactions potentials, used in the literature to model short ranged 
attractions in real solutions of colloids, reverse micelles or globular proteins, can be mapped 
onto this model^. We constructed an approximate excess free energy by allowing various 
coefficients within the BMCSL free energy for hard spheres to be come r-dependent and 
matching to the (for case II, particularly simple) third order virial expansion. The result- 
ing binodal in the monodisperse limit is rather different from the one obtained from the 
PY closure with the energy route. The polydispersity-induced shifts of the (coincident) 
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cloud/shadow curves are nevertheless comparable to those predicted by our PY analysis, 
but only sufficiently far below the critical point. Near the critical point the BMCSL-like 
excess free energy predicts an enlargement of the coexistence region towards higher r, while 
the PY closure gives the opposite result. Given that in the monodisperse case the PY bin- 
odal is rather closer to simulation results than the BMCSL-like one, we would expect that 
also for the polydispersity effects the PY predictions are more accurate. 

The second model for which we considered an alternative approximation scheme was the 
AO model. Here a direct comparison with free volume theory is straightforward since for 
the latter a generalization to polydisperse colloids has recently been derived^. Even though 
one expects the two approaches to be valid in complementary regions (small polymer size ^ 
for the SHS mapping, larger ^ for free volume theory) we found very good qualitative and 
even semi-quantitative agreement of the predictions from the two routes for an intermediate 
value (0.1) of the polymer-to-colloid size ratio. 

In future work, direct simulations of polydisperse SHS mixtures would obviously be of in- 
terest to test our predictions and resolve any differences with other approximation schemes, 
e.g. in case II. Simulations would be ideal here since in contrast to experiment they would 
allow one to probe directly different choices for the stickiness coefficients. Because of the 
presence of polydispersity, a grand canonical Monte Carlo approacb^^j^^i^^i^ may be the sim- 
ulation method of choice, possibly supplemented by specific cluster algorithms tailored to 
sticky interactionsi^i^liiSi. For the physically more realistic AO model, our predictions should 
be more accurate than those of free volume theory for small polymer-to-colloid size ratios. 
Detailed experimental or simulation tests in this regime would be welcome. In simulations 
one could work directly with the AO-depletion potential for the colloids, without ever repre- 
senting the polymers explicitly. For comparison with experiment one would need to work out 
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the actual volume fraction of polymer in the system rather than in a reservoir; this should 
in principle be a straightforward exercise once our excess free energy has been rewritten as 
a function of polymer chemical potential. On the experimental side one would require that 
the colloids are sufficiently polydisperse (beyond a terminal polydispersity around s = 0.07; 
see the discussion and bibliography in Ref.— ) to suppress kinetically any solid phases, thus 
allowing stable observation of the gas-liquid phase splits we have calculated. 

APPENDIX A: PERTURBATIVE EXPANSION OF Lij 

For the perturbative expansion of Eq. (fTHj) one needs the expansions of aij, Pij and (pij. 
These involve the trivial expansions 

(Ti = l + 6i , (Al) 
a,, = l + ^{6, + 6j) , (A2) 
(TiCTj = l + {6i + 6j) + 6i6j . (A3) 

One also needs the expansions to quadratic order of the moments 

Pra = p((l + = p(l + m{5) + ^m(m - 1){5') + ...), (A4) 

giving in particular p2 = p(l + 2(5) + (5^)) and A = l—i] = 1 — ps = Aq — 3p((5) — 3p(5^) with 
Aq = 1 — p as defined in the main text. The final ingredient is the expansion ^ for the e^j, 
which is left in general form to allow different possible choices of the stickiness coefficients 
to be considered together. Altogether one gets the following expansion coefficients for the 



ai„T = (1 + ei„)^ + (I + |ei„) ^ , 

ai6T = 6^ + 9^ , 

Ci2ar = {\ + 2eia + e2a) + (3 + feia + |e2a) ^ 

a;26T = (I + eia + ^26) + (1 + feia + \^2h) -§1 , 

«2cr = (f + 66i„) 4 + (f + 96i„) ^ , 
«2.r = 27^ + f^, 

Similarly one has for the (3ij 

Por = p , 

PlaT = (I + eia) P , 
AbT = , 

P2aT = (eia + e2a) P , 
P2bT = (|eia + £26) p , 
P2cT = , 
I32dr = , 
P2er = , 
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and for the (pij 



^0 - i ' 



- An ' 



1_ 

^0 



k — 3p 
^2e - A? 



One now inserts these expansions into Eq. ()18|) and proceeds as explained in the main 
text to obtain the desired conditions on the expansions coefficients Lq, . . . ,L2e of the Lij. 
To state these, it is helpful to define the quantities 

where Greek indices stand for the labels 0, la, lb, 2a, 2c, 2d, 2e, of the coefficients of the 
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perturbative expansions. The desired conditions are then 

Lo = ao + /5oMo,o , (A5) 

Lla = aia + f3laMo,0 + (3oMo,la , (A6) 

Lib = + Af>Mo,o + 2/3o(Mo,ia + Mo,ib) , (A7) 

L2a = a2a + p2aMo,0 + 2/3iaMo,la + PoMia,la , (A8) 

L2b = a2b + P2bMo,o + f^iaMoM + /5oMo,2fe , (A9) 

L2c = «2c + /?2cMo,o + 2/?i,(Mo,ia + Mo,ib) + AfeMo,!, 

+ l3o{Mia,la + Mia,lb + Mo,2a + Mo,2c) , (AlO) 

L2d = a2d + (32dMo,o + 2f3ib{Mo,ia + Mo,ib) 

+ /?o(2Mo,2c + 2Mo,2d + 2Mia,ib + Mib,ib) , (All) 

L2e = a2e + p2eMo,o + Po{Mia,la + 2Mo,2fe + 2Mo,2e) , (A12) 



The first of these determines Lq and leads back to Baxter's solution ()24|1 for the monodisperse 
case. All other equations involve the desired coefficient on the left at most linearly on the 
right hand side and so are trivial to solve; e.g. Eq. ()A6|) has Li^a on the left and implicitly 
via Mo^ia on the right. Running through the equations in order, all expansion coefficients 
can then be found. 
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LIST OF FIGURES 

Fig. 1 Equation of state, from the energy route, for a one-component fluid of SHS. From left 
to right and top to bottom the four panels refer respectively to a reduced temperature 
of r = 1.00, 0.50, 0.20, and 0.15. The continuous line corresponds to the MSA approx- 
imation, the dotted line to the mMSA approximation, the short dashed line to the CI 
approximation, the long dashed line to the PY approximation, the dot-dashed line to 
the WCA first order perturbation theory, squares to the WCA second order perturba- 
tion theory (with error bars indicating the range where the true value should lie with 
probability 99.7%), and triangles to the MC simulations of Miller and FrenkeF°. In 
all cases the HS component of the pressure was chosen to be the one obtained from 
the compressibility route of the PY approximation^. 

Fig. 2 The overlap volume Vov{r) of the two exclusion zones around colloid particles of di- 
ameter (jj and CTj which cannot be accessed by polymers of diameter ^. 

Fig. 3 Phase diagram of the monodisperse SHS fluid obtained with the PY closure and the 
energy route to thermodynamics. Shown are the binodal and spinodal curves and the 
region where the PY equation has no solution (see Eq. ()25|l ). 

Fig. 4 Pressure from the energy route of the PY approximation for a single (parent) phase 
with case IV stickiness coefficients, plotted against volume fraction. Results are shown 
for several small values of the polydispersity s (see legend) and well above, just above, 
and below (from left to right) the critical point of the monodisperse system. The 
pressure was determined using Eq. (9) of Ref,^-. 

Fig. 5 Cloud and shadow curves for SHS mixtures with polydispersity s = 0.3, as obtained 
within the PY approximation and the energy route to thermodynamics, for coefficients 
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eij chosen according to cases II and IV from Eq. (0). The shifts from the binodal of 
the monodisperse system (labeled "mono") were calculated using Eq. (fT3j) and give 
the leading (0(s^)) corrections in a perturbative treatment of polydispersity. Note the 
collapse of the cloud and shadow curve, as expected to this order of the perturbation 
theory for purely size-polydisperse model o^^i^^ , and the divergence of the perturbation 
theory at the monodisperse critical point. 

Fig. 6 Cloud and shadow curves for the SHS model with polydispersity s = 0.2 and case V 
stickiness coefficients. The binodal of the monodisperse system is shown for compari- 
son. 

Fig. 7 Cloud and shadow curves for the SHS model with polydispersity s = 0.1 and case I 
stickiness coefficients. The binodal of the monodisperse system is shown for compari- 
son. 

Fig. 8 Cloud and shadow curves for the AO model with polymer-to-colloid size ratio ^/o"o = 
0.1 and (colloid) polydispersity s = 0.07. The binodal of the monodisperse system is 
shown for comparison. 

Fig. 9 Fractionation in SHS mixtures with stickiness coefficients chosen according to cases II 
and I, at r = 0.11 and for polydispersities s as in the corresponding Figs. El and [3 
Shown are the cloud (parent) size distribution p{(j), taken to be of Schulz form, and 
the size distributions in the liquid shadow and gas shadow phases that form when 
coexistence is approached from low densities (gas cloud phase) and high densites (liq- 
uid cloud phase), respectively. For case II (main graph) the larger particles tend to 
accumulate in the liquid phase, while for case I (inset) the opposite is true. 
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Fig. 10 Decomposition Aa = Aqq + Cia^c^i of the difference in a between gas and liquid 
phases. The two contributions Aoq and Aai are plotted separately against r; the 
latter quantity is graphed on the vertical rather than the horizontal axis for ease of 
comparison with Figs. El to |H1 Inset: Ratio Aao/Aai. 

Fig. 11 Cloud and shadow curves for case II stickiness coefficients and with polydispersity 
s = 0.3, calculated using the BCMSL-type free energy, Eq. (jSH)), rather than the PY 
approximation as in Fig. The binodal of the monodisperse system, which differs 
from the PY result, is shown for comparison. Main graph: region around the critical 
point. Inset: global view of the results on the same scale as in Fig. El 

Fig. 12 Comparison of predictions for AO model with polymer-to-colloid size ratio $,/crQ = 
0.1. Left: Results of SHS mapping analysed within PY approximation; as in Fig. |H1 
cloud and shadow curves are shown for colloid polydispersity s = 0.07, along with 
the monodisperse binodal for comparison. The vertical axis now shows the polymer 
volume fraction rather than the reduced temperature r. Right: Analogous results 
obtained from free volume theory. Inset right: Fractionation coefficient Aa for the 
two approximation schemes. 
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Tab. 1 Coefficients of tlie perturbative expansion ^ of the adhesion parameters e^j for the 
four cases hsted in Eq. 
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